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Please note the following corrections as you read this Report; 

1. Page 11 - Line after eqn, (15) should read 
Equation (15) is a transcendental equation giving non-harmonic sinewaves 

2. Base 15 - Fourth Itoa balotr tha aquation , raplaoa aatandei by 

attended* 

3. Page 21 - First line below Eqn. (42) should read 

and the elements tr°^ are also uncorrelated ... 

4. Page 21 “ Eqn. (42) should be 

~ ^1®2 ^ij^ikSj^ (^2) 

5. Page 23 - Fifth line, second paragraph should read 
..... subtract the image mean before ... 

Second last setence should read 
Thus, once again x.e see that even though the covariance function of u. . 
is nonseparable, the ICL transform of given boundary conditions 
is a fast transform, i.e., the sine transform. 

7. Page 26 - Fourth line from bottom should read 

^i ~ ‘ the i component and ... 

8. Page 28 - Ninth line from top should read 

Thu bnundury Vg. , -I 

9. Page 28^- The equation for near the bottom of the pages the firat line 
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a* f 2 
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10, Page 30 - Fourth line from bottom should read 

scheme is . . . . with g 0.7 , this was . , , . 

11. Page 23 - Line 15 . * , . bock- vblock 
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ABSTRA.CT 


The purpose of this study was to apply the fast KL transform algorithm for 
data compression of a set of four ERTS multispectral images and compare its 
performance with other techniques studied by TRW, Inc. on the same image data 
in contract NAS2-8394. The performance criteria used here are mean square 
error and signal to noise ratio as used in the TRW effort. The results obtained 
here show a superior performance of the fast KL transform coding algorithm on 
the data set used and with respect to the above stated performance criteria, 

A summary of the results is given in Chapter I and details of comparisons and 
discussion on conclusions are given in Chapter IV. 
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CHAPTER I ORIGINAL PAG-B IS POOR 

INTRODUCTION 

1. Image Transforms 

Two dimensional image transforms have attracted considerable attention 
during the past decade for their application in image data compression (commonly 
called 'transform coding'). If U represents anN x N image matrix, then, in the 
context of image transform coding, W is its two dimensional transform where 
W = TUT' (1) 

where prime denotes the transpose and T is a unitary two dimensional N x N 
matrix. The restriction of T to unitary matrices ensures conseryation of image 
energy in the transform domain. 

Examples of image transforms are Fourier, Hadamard, Karhunen Loeve (KL) , 
Haar, Slant, Cosine, etc. Typically, the transformed image W is such 

that most pf the image energy is concentrated in relatively few samples in the 
W-space (usually the lower 'frequency' samples) so that these few samples alone 
are considered important for any subsequent image processing thereby achieving 
some data compression, 

2. Dimensionality vs Optimality 

Tvjo considerations which become important in selecting the image transform 
for data compression are a) Dimensxonality and b) Optimality for Compressxon. 

By dimensionality we refer to the computational effort required in implementing 
Eqn. (1). For an arbitrary unitary transform T, it is easy to see that compu- 

•- *3 

tation of W takes 2N multiplications and about as many additions. For images 
the value of K gets large enough (N ci 200 to 1000) to make this computational 
load unbearable. Consequently, the evil forces of dimensionality dictate that 
the choice of T be restircted to certain class of transforms which can yield 
fast computational algorithms* One class of such transforms called Good Trans- 


rx 81 

forms (after I. J. Good^ have the property that the matrix T can be written 
as a product of several sparse matrices in the form 

T = TiT2...Tp (2) 

where T^, i = 1, , , . ,p Cp<<N) are matrices with just a few non-zero entries 

r entries per row, say, with r « N) . Thus the multiplication of T with, 

say, an N X 1 vector, is accomplished in about rpN computations. Thus, if 

N =! 2^, then typically, the computations required in (1) reduce to approximately 
2 

rN log^N. Depending on the actual transform, one computation can be defined 
as one multiplication and one addition/ subtraction, e.g. , Fourier, Cosine trans- 
forms or as one addition/subtraction, e.g. , Walsh-Hadamard transform. 

By oftimality we mean the efficiency of a transform in achieving data 

compression (or bandwidth compression) . Usually this optimality is measured 

% 

for a class of images rather than for a single image because, conceivably a 
transform could be optimal for one single image and be very poor for others. 

This raises the question, "How do we define a class of images?" Here we con- 
sider classification of images by their statistical properties as opposed to 
their linguistic or descriptive properties, although many times it may be pos- 
sible to quantitatively model linquistic/descriptive images by a statistical 

expression. For a given class of images having certain second order statis- 

’ r 9-13 1 

tics, the Karhunen Loeve Transform* is shom to be optimal in a mean 

square sense. 

Although the KL transform is optimal, it has dimensionality difficulties. 
First, the KL transform is unique for a class of images. Therefore, it has to 
be computed for that class. Second, even if a closed form analytic expression 
for the KL transform is known, the transformation calculations of Eqn. (1) do 
not, in general, have a fast algorithm available. Therefore, the size of image 

F7 10 1 

* Sometimes also called Hotelling Transform after Hotelling.'- ’ 


tha-t can be used for KL transform application is quite small 8 x 8 or less) , 
An example of utilizing the KL transform technique is by dividing an image into 
small blocks and then coding each block. Otherwise, most transform coding 
efforts have resorted to the class of fast image transforms described above. 

3. Discrete Cosine, KL, and Fast KL Transforms 

Recent experimental results of Ahmed, Natarajan and Rao [5] have shown the 
behavior of the Discrete Cosine Transform (DCT) to be close to the KL trans- 
form for one diraensioanl first-order markov random porcesses whose correla- 
tion parameter p has a value around 0.9. Actual implementation of the DCT 
for 2-dimensional image coding has shown that the DCT does perform better than 
Fourier, Walsh-Hadamard transforms, etc., at low bit rates (bit rates of less 
than equal to 1 bit/pixel average) . However, at higher bit rates (n^ 2 bits 
per pixel average) , a hybrid coding algorithm employing the DCT in one dimension 
and DPCM in the other seems to perform better than the 2-diraensional DCT |*20]. 
This, of course, should not happen if the DCT were truly a good approximation 
to the KL transform, because tHe rate distortion curve of the KL transform is 
expected to give a lower bit rate at a given distortion level. This ambig- 
uity regarding the 2-dimens ional DCT may be partly be attributed to the fact 
that image statistics is not gaussian, and partly that the image covariance 
function is not separable, 

4. Summary of Results 

In this study, the fast KL transform algorithm developed in the next two 
chapters is applied to a set of four channel ERTS multispectral images and 
its performance is compared with other techniques studied by TRW, Inc. [20] 
in contract wAS2-8394, The image data used in this study is identical with 
that used by TRW, Inc. and in fact was supplied to us by them* Figure 1 


summarizes the rfesults of this study. Here each curve represents average 
of mean square error over the four channels vs. fait fate obtained by a given 
technique. Thus if 

= mean square error in the k channel image, k— I 3 • . .4 (2) 

then 

^2 

average mean square error - % 2 e-i^ 

k=l 

The KL-Cosine/DPCM, I<L-Hadamard/DP(M curves, correspond to three dimen- 
sional data compression schemes utilizing a 4x4 KL transform along the temporal 
axis. A.11 curves in Fig. 1 except the 2 -dimensional fast KL curves were obtained in 
the TRW study C^O]. The fast KL curves correspond to two different sets of 
experiments; viz, the first set (solid line) using 255 x 255 size transfor- 
mation and the second set (dotted lines) using 15 x 15 block transformations 
on the image. In the other curves 16 x 16 spatial block size was used in per- 
forming any transform coding. Comparisons show a superior performance of 
the fast KL transform coding algorithm over other methods. 

A detailed discussion of these and other comparisons with respect to 
signal to noise ratio etc., as well as conclusions of the study are given in 
Chapter IV. 


Bit Rate - Bits/ Pixel 



FIGURE 1 Coiaparison of mean square error for 2 Dimensional Fast KL vs. Other Methods 
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CHAPTER II 


A THEORY OF FAST KARHUNEN LOEVE TRANSFORM ' 


1. Image Covariance Model 

An image may be tbought of as a sample function of a two dimensional 
random process. If u. . denotes the brightness at the spatial coordinate (i,j) 
then u is considered as a random process . Consider one such process which 
has zero mean unit variance and is described by a covariance function 
given by 


RCnjm) = E[u 


i+n, j+m'^ij ^ 





where 


1 > p .= correlation parameter j i - 1, 2. (5) 

The covariance derihed above assumes essentially two properties viz., 
stationarity and separability. Although one might aigue for an isotropic non- 
separable covariance function of the form 

2 ji?-hn^ (&') 

R(n,m) = o p«l lo; 

as being more appropriate, the covariance of Eqn. (4) has been found to 
be a reasonable assumption and has bee:n applied successfully in many different 
image processing problems . ^ ^ The primary reason for preferrijig Eqn. (4) 

over Eqn, C6) is due to the amenability of covariance in (4) to simpler 
analysis and applicability of many one dimensional results. 

2. One Dimensional Representation 

Let {x.} be a finite one dimensional random process with zero mean and 
X ' 

unit variance and a covariance function given by 


8 


E[x.x., ] = i = 0,1,.. ,E+1 . (7) 

1 i+n 

The zero mean and unity variance assumptions are non-essential hut serve only 

to present a simplified analysis. It is well known that the sequence x* can be 

represented by a first order stationary Markov process as 

X. ,, - px, + e. i - 0 (8) 

r+1 X X 

with ■ 

E[e.] = 0 

X 

and 


E[e.e,] = (l“p )5.. , 6.. = Kronecker delta function. 

X 3 xj ’ xj 

It has been shorn earlier Xl, 14] that the above Markov process for a fixed 
N, can also be represented by . . 


X. = a(x.._+x. _) + V. 
X x-hl x-1 X 


1 ^ i ^ N 


( 10 ) 


where 


or - 


1 + 


( 11 ) 


and i = Ij • • »3 Ef} is a well defined random process. 

The above representation is defined for a finite process (N = fixed). 
The interpretation of this representation is as follows. If x^ denotes the 
best linear mean square estimate of x^ obtained from a linear combination of 
the elements of the partial sequence [x., jj-i}, then the residuals {v-} 
given by 


X. - = V. 

X i X 


9 


are such that the variance 


E[(x^-x^)^3 = E[v^^] 

is minimum. This will be called the minimum variance representation. In 
order to find we start by writing 


N+l-i 

X. = I a., X. , + E b., X. 
1 lt=3 |c=l 


i = 1, . . . , N 


_ 2 

and find a., and b., for each i and k such that E[(x.-x.) ] is minimized, 
xk xk XX 

l^hen this minimization is pei'formed [14], the coefficients a^^ and b^j^ are 
obtained as 


a. , 

“ a - b , , 

VI 

H 

xl 

il 


a . 

= 0 = b 

k - 

xk 

ik 


The 

’v' Process 



k-2, 1-i^N 


The minimum variance property of the sequence does not guarantee its 

being uncorrelated. In fact, in our case, the sequence elements hrxve nearest 
neighbor correlation. The corralation properties of the sequence can be 

expressed as 


E[v.v 1 

X J 


1 ° ■ 

for i,j = 1 N, and 




i = j 


li“il = 1 , P = 2 

1 + p 

otherwise 


(12a) 


Erv^x.] 


(12b) 


These results are derived in Appendix I. 
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4. The Karhunen Loeve Transform of Markovian Images 

If X is a one dimensional n x 1 vector with covariance matrix R, 
then the KL transform of x is a matrix composed of the eigenvectors of R 
and is defined by the relation 

$’Rf == r (13) 

2 

where r is a diagonal matrix of eigenvalues If x if a first-order markov 

process with covariance given by (see Eqn. (7)) , then 
, 2 n-1 

1 p p p 

^ n-2 I 

plpv^-- — p 


p^'- p -1 , 


\ 


R = 


i , • ^ 

i , p 

i n-1 2 , 

;P P p 1 


and the elements of the KL transform are given by [11] (for n=even) 


, . r / » n*b 1\ , 

. = a. sin [oj. (j TT-J + -r-J 

ij 1 1 2 2 


(14) 


where y . = 
' 1 


1 - p' 


1 - 2p cos tii^ + p 


77 , a. is the normalization constant and {m,} are 
2 ’ i 1 


the positive roots of 


tan noi = - 


( 1-0 ) sin m 


(15) 


cos (i) - 2p + p cos (jj 

Equation (15) is a transcendental equation given non-harmonic sinewaves 
of Eqn. (14). The KL transform of the vector x may be written as x = 4>'x or 


n 


n 


x.= Z i)..x.= I a.x. sin [tD,(i + ) + '^] 

^ j=i J 1-1 ^ ^ J ^ ^ 


(l6) 


Due to non-harmonic behavior of the sine terms, a fast algorithm (like the EFT) 

2 

is unavailable in computing the series of Eqn, (16). Therefore, typically n 
computations are required in computing i “ l,...,n). 


U 


RijPEODUcrBiLrrY of the 
Smal page is poor 

From (13) and (16) and using ©' f = I, it is easy to see that 

. - , ,2, (17) 

Et;x,x.] = Yi6. . 

Since the samples are uncorrelated, they are quantized independently. 

Other advantages o£ the KL transform are its minimum entropy and minimum mean 
square error properties which make it optimal for data compression [9, 13]. 

These properties assure that for a chosen ratio of compression, minimum mean 
square distortion will result; this in comparison with all other linear uni- 
tary transformations for the assumed class of signals. 

5. The Fast KL Transform for Markovian Images 

Consider the sequence k - 1,...,N], and represent this by a vector 

V, then the N x N covariance matrix C given by (12) is 


-a 


-a 


-a 


C = 3[vv’ 1 = 3 


-a 


(18a) 


-a ' 1 


i 8^ Q 


(18b) 


tri- 


where Q is the N x N matrix in (18a). The matrix Q is a symmetrxc, 
diagonal, Toeplits matrix. 

Theorem 1: The ICL transform of the sequence (vj^, k = 1,...K] is given by 


2 . ill 

'I'ij =Jn+I "^^n+1 


( 19 ) 
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Proof : 


The eigenvectors i]). . and the eigenvalues X , of the N x N symmetric, tri 
diagonal, Toeplitz Q are given by [14] 

4 




2 . 1,1 TT 

ij N+1 N+1 


( 20 ) 


and 


1 T o 3.H 

X. = 1 - 2a cos ^ 


for i, j = 1, . . . ,N 


( 21 ) 


Clearly, the matrix C in Eqn.(18a) has its eigenvectors also given by 
(since 3^ is a scalar constant), so that the matrix is the KL transform 

of V. 


Theorem 2: [1] For the first order, stationary, finite, Markov sequence 

{x^, i = 0, 1,...,N, N+1}, if the boundary conditions x^ and are given, 

then the KL transform of the partial sequence {x^, k = 1,...,N} given x^, x^^^^ 
is the matrix * vrith elements given by Eqn. (19) . 

J-J 


Proof: Let the given boundary conditions be 


^0 ~ ^+1 


= d 


( 22 ) 


If X and V are defined as N x 1 vectors of components {x^,...,Xjj} and 
respectively, then Eqn. (10) can be written as 


Qx = V + b 

where Q is the N x N tridiagonal matrix in (18a) and b is an N x 1 vector con- 
taining only the information at the end points; viz . , 


b, = ac, b„ = ad, b, = 0, 2 ^ k ^ N-1 


(24) 


^1 “N k 

Since ’c’ and ’d’ are now given, and v and b are uncorrelated (see (12b)) 

4 _ , I .I ^*"1,^ ( 25 ) 




^ y = E[x|c,d] = Q H 


and 


T,^-l _ „2„-l 


= E[ (x-y) (x-y) lcj,d]=Q E[vv]Q =8Q 


(26) 
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2„-l 


Since 


Hence, the covariance of x given end conditions is simply — 3 Q 

2 2 "1 
S is a scalar, the eigenvectors of p Q are given by defined above. 


Equation (23) above can be rewritten as 
x = Q v + Q D = x + y = x 


(27 ) 


where we have defined x = Q v and Xj^ is defined in (25). Thus, Eqn. (27) 
decomposes the random process x into two terms, viz., x° and x^ where x° has 
zero mean and covariance given by Eqn. (26) and Xj^ is the conditional 
mean of x given boundary conditions. Denoting x = ijjx and similarly Vj S, x°, etc. 
and realizing from theorem 1 that = A, Eqn. (27) can be transformed to 


yield 


V. b. 

A 1 , 1 A O , A 

X. = — + ” = x. + y. 
1 X. \ . X X 
X X 


(28) 


The definition of b in (24) gives 


A ( 2 

\ =>/ 


N+1 


a(c-(-l)^d) sin 


ITT 

N+1 




(29) 


and application of the result in (26) shows fx?'! are uncorrelated, i.e.. 


-AOaO.| jK 

E[x.x. ] - 6 


IJ 


(30) 


From Eqn. (26), the eigenvectors are independent of the correlation para- 
meter p and only the eigenvalue depend on the statistics of the random pro- 
cess x. This is in contrast with the eigenf actions 0.., Eqn. (14), which depend 

2 

on p through and Moreover, the eigenvectors of Eqn. (20) are harmonic 

sine waves, so that a fast computational algorithm is possible and is developed 
in the next chapter. 

The decomposition in (27) has some interesting properties. First x, , the 

b 

conditional mean of "x given boundary conditions," depends only on the boundary 
conditions, as given by (25,). The process x^ is called the boundary 


14 


response of the process x. The process x° is obtained aiaiply by subtracting 
the boundary response from the original random process x. Second, the decom- 
position is orthogonal, i.e., 

E[x°xJ] = 0 (31) 

This follows from Eqn. (5) of Appendix I where 

E[ViXi^^] =0 0, l^i^N 

Since x° = Q ^v, and x^ contains Xq and only, E[v^Xq1 = 

When the boundary conditions are known (i.e,, c, d are given) then Xj^ can 
be easily calculated and the vector x° is written as a modification of the 
original vector x according to 
o A 

X = X - x^^ 

The vector x° can now be coded by its KL transform, which is a fast trans- 
form. The original vector x is recovered easily from x° if the two boundary 
values c, d are separately coded for transmission. This idea above easily 
estended to two dimensions when the mo dimensional image autocorrelation is 
separable. 

The special properties of the fast KL transform are: 

(i) The fast KL transform is independent of the image correlation para- 
meter, p. Only the boundary response x^ depends on p. Hence, the fast KLT 
coding algorithm could be useful in adaptive coding schemes, compared to the 
actual KL transform which varies with p. 

(ii) The transform domain variances are known as a closed form expression 
for any fixed p. Hence, the quantizer design calculations are simplified (since 
the transform domain bit assignments depend on the distibution of these vari- 
ances) . These calculations are specially facilitated if the quantizer is to 

be adaptively changed with the changes in correlation parameter. 
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Ciii) The number of computations required in fast KT transform calcula- 
tions for an N x N image is of the order N log 2 N> the same order as in FfT, or 
fast DCT. Strictly speaking, the fast KL transform computations are less than 
the fast DCT computations. 

It should be noted that the fast KL transform is not a numerical algorithmic 
fast solution of the conventional KLT of the data, rather, it is the conventional 
and fast KL transform of a modification of the data. Figure 2a shows how this 
modification is achieved. The boundary values , and are passed through 

a linear filter to compute the boundary response Xj^(l) , . . ., Note 

that this entire response is generated by only two input values. The modified 
data is then x°(k) , 1 £ k ^ K, which can be fast KL transform coded. The 
boundary conditions shown to be uncorrelated with x (k) , and are 

therefore coded independently. Figure lb shows the reconstruction of the data 
from the fast KL coded data and the boundary conditions. The linear filters 
employed in figures 2(a) and 2(b) are known apriori and each filter requires 
only 2N computations for a N x 1 vector for computing its output. The filter 
calculations and coding of boundary values can be avoided in some practical ap- 
plications by approximating the boundary values by the mean value of the data 
vector X and subtracting this mean value from the data vector before compres- 
sing it via the fast KL transform. Hence, for zero mean data the boundary con- 
ditions may be approximated by zero and only the mean value of the data vector 
is coded in addition to the modified vector x°. 


6. Two Dimensional Representation 

Let the image be described by a zero mean random process with i,j = 

0,1, ...,N,N+1 and autocorrelation of Eqn. (4). Following Appendix II it can be shown 
that the minimum variance representation of u^^ can be written as 


u, . 


a(u. 


.+u . 


i+l,j “i-l,j 




V. . 
11 


V. . 
11 


- a(v. ..,+v, . ,) 

2 i,j+l 1,1-1 


V. . 

ij 


B^RODUCmmOT OF 'iffl 
original PAGL is iOOR 


(32) 

(33) 
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where 


•p-GjppQX)UGIBILTrY OF 

oSSl page is pgg» 


i.J - 1,2,. 


, N and cs?j^ = 




k =1, 2. 


If v_ is eliminated from Eqns, (32) and (33), we get 




■** 1 .5+1 1 J 1 ) + V . , 

i-J-jJ+l x-1,3— 1 X3 


(34a) 


- "ij ^ij (34b) 

Figure (3a) shows the as a linear combination of the nearest eight neighbors 
in this case. 

7 • The Two Dimensional Fast KL Transform 

The two dimensional result is summarized by the following theorem. 

Theorem 3 : If u^_. is a finite, tiro dimensional random process defined for 

i,j = 0,1, . . . ,N,N+1, with its autocorrelation given by Eqn. (4), and if the 
boundary conditions for all i,jE[0,N+l]} are 

given, then the KL transform of the partial field (u , k = 1,.,,,N; = 1,...,N} 

is given by 

T = , 

where ijj is SINE transform defined by elements 


'^ij i+1 

For proof of this theorem we reisort to the two dimensional minimum variance 


representation of u^^ in (32) and (33) . If U,V, and v are N x N matrices of ele- 
ments u.., V.., and v . respectively, then 
aj ij ij 


;i8 





Q^, U == V -5- 

VQ^ = V + 
where and 
are given by 


( 35 ) 


are N x N matrices containing only the boundary information and 


h 


^ 0,1 “ 0,2 


u, 


“n+1,1 


0,N 


<— 1st row 


N-2rows of 'zeros* (36a) 


I ^ 

^N+-1,N|<-Kth row 


1st 

column 


V 


1,0 


V 


2,0 


®2 


Nth 

column 

^1,N+1 ‘ 


2,N+1‘ 


(36b) 


!"n,o 


N-2 columns of zeros 


, "^N,N+1 
-)1 — 


In Eqn. (36b) , the quantities v. and v. are obtained from the boundary 

1 3 O X y iM'J’X 

information of u. and u. for i = 1,...,N, via Eqn. (32). Eliminating 

x,o ijN+1 

V in Eqns. (35) we get 

QjUqj. V -i- 

= V + B 

where 


( 37 ) 


B ^ B^^Q + B^ 
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contains only the boundary information. If U, v, and B are now defined as 
2 

N X 1 vectors corresponding to a lexicographic (dictionary) ordering of the 

T 

matrices U, \i and B respectively ^ Eqn. (37) can be rewritten, after noting as 

(Q^© Q2)U = v + B (38) 

where © denotes Kronecker product. From Appendix II, the vector v has zero 

mean and E[ijv'^] = Hence, from Theorem 1 in the last section, (i{i©i|)) 

2 2 — 
is the N X N KL transform matrix of v and 


- p = E[u|b] = (Qj®Q^)"^ 


(39) 


and 


= E[(U-y) (U-iI)^lB] = P^32^Q^^Q2)'\ 8^ = 1 , k == 


- 2.2 


-1 


1-p, 


k 


1+Pi 


(40) 


Hence the matrix of eigenvectors of also is (\jf © iji) , so that the KL 
transform of given boundary conditions h, £ = 0, N+1} is given 

by the matrix ij; defined in Eqn. (20) . Finally, if Eqn. (37) is pre- and post- 
multiplied by t|) and we define W = iJiUij/ , v = , B = ijjBijj , the following rela- 
tions are easily verified, if X^(i) eigenvalues of Qj^, Q 2 respectively. 


w. . 

13 


V.. AO 

# - “ij + hi - h<-» 


- 4 44 «ik 

and the elements U.. are aslo uncorrelated. 

• • 2.J 

orthogonal decompositon of the image matrix as 


(41) 

(42) 


Again, from (38) we get the 


U = U° 4. U 


where ti = (Q^ © " '^^1 ® 

or in matrix form 

U = U° + U. 

b 

where U° = q“^ V = q"^ B q"^ 


(43a) 

(43b) 

(44a) 

(44b) 
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(45) 


P,E^TtOBiJCIEHJa:Y C? TH3 
r,-- 'kV' ''-li! H5 


If we define 


it 

X. . 
3-J 


= 2!ii 


A A iV 0 0 

then E[x„] = Oj E[^x.^Xj^^] Thus are identically distributed 

A 

variables of equal variances. If x. . are quantized, then the dynamic range of 

A. 

each quantizers (corresponding to x. .) will be the same although the number of 

xj 

bits assigned to x. , will be different. 

XJ 

8 . Comments : 


The fast KL transfrom has two steps. First, the residuals v- • are calculated 

^ J 

and then a sine transform of the residuals ' is taken. Alternately, (see 41) , 

A 

we can first take the sine transform to obtain w, . and then subtract B. ./x. . = 


3-J 


XJ XJ 


(Uj^(i, j) = the sine transform of the boundary response) to obtain U — . In 

A A 

any event, the variables (or scaled variables x^^) quantized and coded in 

the sine transform domain. In conventional sine transform coding, on the other 

hand, the sine transform samples w. . will be coded. Further, the variances of 
''O 01^2 

U . . are which are different from the variances of U. . in the sine transform 

J.J X. . 

XJ 

domain. Thus, the quantizers in the sine transform domain and fast KL trans- 
from domain will be different. 

Comparison with Discrete Cosine Transform (DCT) 

The DCT of a sequence is defined as 


/2 2 (2nri-l)kTT 

= E X,g, = “ E X cos 

N m’ ^k N m 

m=o m=o 


2N 


, 1 i k - N - 1 


(46) 
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A-lthougli both the fast KL and the Discrete Cosine transforms have obvious 
relationship with DFT, the DCT does not satisfy the KL transform Eqn. (13) for 
either R (the original covariance of Eqn. (4) or for R^^ (the conditional auto- 
correlation of Eqn. (26)). The relationship between the KL transform of Eqn. 

(14) and the fast KL transform of Eqn. (20) is more easily seen; the eigen- 

xtt 

vectors in the latter being periodic with of (15) being replaced by . 


Necessity of zero mean image 

It is easy to check that the sine transform does not diagonalize a 
matrix of all constant elements. Thus, if the sine transform 
is the KL transform of a sequence jx°, 1 s i ^ N}, it is not the KL transform 
of the sequence £x^ + a^ where *a* is a constant. Since the given image data 
may not have zero mean, it is essential to subtract the image mena before ap- 
plying the fast KL image coding algorithm. In blc_ block coding it is 
advisable to make each image block zero mean before coding. The overhead of 
transmitting the bock mean for 15 x 15 blocks, 7 bit image is only 0.031 bits/ 
pixel. 

g . Ext ension to Markov Images with Nonseparable Covariances Eunetions 
A random field described by the stochastic difference equation. 


u.. = (u. - . + u . . + u- 

XJ 2; J 3 XjJ X Xj J*ri- Xj 

where E[e.^>0, 

-0? k = +1, f = 0 
2 




0 Otherwise 
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(47) 

(48) 

(49) 


generates a Marlcov-1 field [21] whose covariance function is nonseparable. 
Such fields could give a better approximation to the actual covariance func- 
tion of images. If an image is represented by an N x N segment of this field, 
then following the notation of the previous sections, we can write 


QU + UQ = e + B (50) 

where Q, U, B, e are N x N matrices; the matrix Q being tridiagonal as before, 
and B containing only the boundary information (this B matrix is obviously 
different from that in (37)). In Kronecker product notation, it can be seen 
that 


E[ei^] = (I © Q + Q 0 I) = R 
A - „ 

B = E[U 1 B] = R B 

R^ = E[(U -u ) (U -rt )^lB] = R 


(51a) 

(51b) 

(51c) 


rh i = j 

where Q..=j - | |i - jj =: 1 
' 0 Otherwise 

The KL transform of R is i|t 0 \ji aud hence \{i 0 ijt is also the KL transform of 
R^ in (51c). Thus, once again we see that even though the covariance fucntion 
of given boundary conditions B, is a fast transform, i.e. the sine trans- 

form. Hence, images represented by random fields of (47) can be fast KL coded 
by calculating p, (boundary response) , subtracting it from U and coding the re- 
sidual U - p via the (fast) sine transform. 
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CHAPTER III 


THE fast KARHUNEN LOEVE TRANSFORM CODING ALGORITHM 


1. Implementation on One Dimensional Image Data* 

It should be pointed out that while Eqn. (20) is directly related to the 
FFT, the elements are not the ’imaginary' term in the conventional Discrete 
Fourier transform. It is easy to show that ij; _ form a complete orthonormal set 
of basis vectors, whereas the imaginary terms of a DFT do not form a complete 
set of basis vectors. 

Let 

A N r-r^ K 

Jl 


y = £ ’f'l ^ sin 


'k ,:i^k££ 


N+l 


)l=l 


N+1 


If the DFT is defined as 


M-1 


2nik 




0 < h s M - 1 


(53) 


then, in order to implement (52) define a variable as 


z = 0 
o 




z^ = 0 


, i. = 1,...,N 

, £ = N+l,N+2, .. ,,2N+1 


(54) 


Now let z, be given by 


A 

Z,. = 


2N+1 


^ i“2(N+l) ®-=0 


£ z„sinfg = Im[DFT{z^}] = ^ 


(55) 


Sinte z can be computed via FFT, y. is obtained via this fast algorithm. 
<. k ^ 


* The x,y, z variables in this section are not the same as used elsewhere. This 
section only shows the algorithmic implementation of the sine transform. 
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Data Compression Via the Fast Karhunen loeve Transform 
First consider the one dimensional case. From (28) we have 


V- b . 

\ ^i 


b. 


where, b^ is known whenever c and d are known. The next question is 'What 

A A 

is the bit assignment strategy for quantisation?' From Eqn. (56), if z.and b. 

A 

are observables, the original signal x can be constructed by calculating x., 
where 


" _ '^i 

"i" Ai ’ 


] = b'^6. . 


Let c*, d*, z* = quantized value of c,d, z. respectively 

X 2. 

^ A A 

xy = reconstructed value of x. from z*, c*, and d*. 


The expected mean square quantization error can be written as 


A . 

1 

(_l)i(d_d*)]^] 

Letting c* = c and d* = d, we get 

A 

1 q. 

‘i 1 = h w-i- ] - r 


If q^ is small, then to minimize Z the number of quantization levels for 

^ 1 

the ith component z. should be proportional to — . 

^ ^i 

Let n. = number of bits to be assigned to the ith compenent and let 


Oi = bi - b2 log2X^ 


if p' = desired bit rate in bits/pixel, then 


p’N = Z n^ = Nbj^ - b2 r iog2A^ 
i=l i=l 
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ITT 

Now recall that X. = 1 - 2a cos rrrvs 
X N+1 


Clearly, 0 ^2 ^N ^ for all 0 < p < 1, 

Hence, • . . > 

Let “ maximum number of bits to be assigned be fixed; then 

^1 - bi - bj logjij 

A.fter solving for b^ and b^ we get 

A, 


n. = n, - 
X 1 


(n^-p')N log ^ 

N Z. 

S log ” 
k=l \ 


(58a) 


The result of the above equation will not give integer values to n^, in 
general. So a truncation to the nearest integer has to be made. 

Let n* = [n,] = nearest integer to n^, then the actual bit rate is given 

by 

1 ^ 

^ N . 1 1 

1=1 

Proceeding as in the above one dimensional case and assuming negligible 
quantization error in the boundary information and defining 
n* . = number of bits to be assigned to v. , 

p' = desired bit rate in terms of bits/pixel 


we get 


n*. = [n. . ] 
11 iJ 


(n --p') log 

.,2 JOI 

^11 “ ^ N N 


(58b) 


hi 

j Z Z log — **- 

i 1=1 j=l ^11 

I- -j 

where [ ] indicates the nearest integer. The two dimensional coding algorithm 
can be implemented in the following steps. 
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3. Fast Karhunen Loeve Transform Coding Algorltlm 

1. Calculate the statistics (i.e.. mean, variance, horizontal and vertical 
correlation parameters) of the source image. 

2. Create zero mean image, by subtracting the image mean (or mean of the 
image block in case of block by block coding) from each point in the image. 

3. Use Eqns. (34a) and (34b) to compute u^^ and for i, j = 1,...N. 

4. Take the two dimensional sine tr'^nsform of v. . to obtain v... Calculate x. . 

xj xj xj 

5. Quantize using n*^ bits to obtain x*^ . Details of quantization 
method iiX& given in the next chapter. The values x_ are used for transmission. 

The boundary information, i.e., {u. , u, „ . , , u u ,, i, j = 0,...,N+ll 

is quantized separately. In practical situations this information may be 
assumed to be the mean value of the image block. 

6 . At receiver compute B*. from {u* , u* „,,, u* u*,, If r , 

xj i,o’ i,N+l’ o,j N+ 1,3 0,3 

defined as the sine transform of the row {u* i = 1 ,...,N}, i.e. 

0 , 3 ’ -J > > 


is 


A 

r 


o »3 


^ r» 4 . *3 K.7T 

N+1 N+1 ’ 


A 



3 = 1, . . . ,N 


A ^ 

and r ., c. , and c. are similarly defined as the sine transforms of 
N+i ,3 1,0 i,N+l •' 

Iu*^l j^}, {u* ^}, and {u* respectively, for i, 3 , k = 1,...,N then it 

A 

can be seen that the boundary terms B*. are obtained as 

X3 

S ('cj + iw 


N +1 '03 ' ' "N+1,3' N +1 

3 + 1 a ri ' . XTT . .iff 

/i5T 

V 1+1 .. . / ,%i+l A j. 


(c +(-i)j+^a. 

xo i,N+l 


[u*o + "^+1,0 + ( ^ '^+l,N+l 


] 


X* . . 

7. Compute wi?. = ^ and obtain u*., the reconstructed image 

/Aij X 3 


X3 


at the receiver as 


[u..] - U* = 

Figure 3b shows the block diagram flow of the above algorithm. 
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4. Discussion 


In practical implementation of this algorithm, N + 1 should equal 2^^ for some 
integer n, so that an FFT algorithm can be used to implement the sine transform, 

The value ofn^^ in (58b) is an experimental parameter to be chosen to minimize 
the mean square error. In our experiments we found the values of n between 6 
to 8 quite satisfactory. These values were determined by testing the mean square 
error of a few randomly chosen image blocks of Channel 1 image. 

For large image block size N ^ 63, and/or small correlation parameter value 
p s 0.7, the effect of boundary values on coder performance becomes less signi- 
ficant. Hence, instead of taking an (N4-2) x (N+2) size image block and using the 
actual boundary values , if an (NxN) image block is used and all the boundary values 
are assumed to be equal to the block mean then the performance of this coding 
scheme is unchanged. For 15 x 15 image blocks with p s .07, this was experimentally 
verified. Note that this approximation does not change the actual boundary value 
around each image block, but only uses the mean value as the pseudo boundary of the 
blo'ck. 
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CHAPTER IV 


FAST KARHUNEN LOEVE TRANSFO^ DATA COMPRESSION STUDIES 

The algorithm outlined in the previous chapter was simulated on four 255 
X 255 picture elements multispectral earth data images. The digitized images 
were provided by NASA and are the same as used by TRW, Inc. in a 
data compression study under Ames Research Center contract NAS2-8394. m all 
coding experiments, the boundary conditions were assumed to be constant and equal 
to the mean value of the image block being coded. 


1. Source Data 

a. Figures (4) through (7) show the original four multispectral earth data 
images. 

b. Table 1 lists the calculated statistics of the original four multispec- 
tral earth data images. 

c. Figures (8) and (9) show the histograms of the original Channel 1 and 
Channel 3 multispectral images. The histograms of Channel's 2 and 4 
are similar to those of Channels 1 and 3 respectively. 

In image coding simulations, two different block sizes were used. First, 
the entire source image was used as a single 255 x 255 image block. Next, the 
255 X 255 source image was divided into 289 (15 x 15 ) small image blocks. Un- 
less otherwise stated, (see Tables) uniform quantizers were used in all the 
coding experiments. The improvement of a compandor over a uniform quantizer was 
found to be less than 0.5db in terms of S/N ratio. The performance of different 
coding schemes is judged by (i) Mean square error and (ii) Signal to Noise Ratio 
This is to enable us to make comparison with other techniques. The mean square 
error in encoding the kth channel image is defined as 


= i 

n2 


N 

S 


N 

E 


i=l j=l 




i j j } k 


1 *. . ,) 


k = 1,2,.. .,4 
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FinURF. 5 CHANNEL 2 I 


FIGURE 4 CHANNEL 1 IMAGE 


FIGURE 7 LHANNE 


FIGURE 6 CHANNEL 3 IMAGE 




Statistical! 

Parameters ^ 
1 

i 

1 Channel No. ; 

: 

Mean 

\ 

i 

Variance 

1 

Average Correlation Parameter p 
(255 X 255 Coding) 

Average Correlation Parameter p 
(Block by Block Coding) 

p hori p ver 

channel 1 

47.80 

• 130.76 

0.93 

0.710 

0.597 

Channel 2 ; 

52.13 

333,16 

0.92 

; 0.730 

0.620 

Ci.anne 1 3 

; 

59,56 

121.16 

0.88 

' 0.712 

0.591 

1 

Channel 4 ' 



26.82 

41.10 

0.89 

0.723 

0.614 


TABLE 1 The Statistics of the 4 Multispectral Earth Data Images 
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FIGURE 8 HISTOGRAM OF CHAITNEL 1 IMAGE 
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FIGURE 9 HISTOGRAM OF CHANNEL 3 IMAGE 
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where N = 255, u. . , and u*. . , represent the original and the reconstructed 
i,j,K i,j,k 

images for the kth channel. Channels 1,2,3 are given to be 7 bit images and Chan- 
nel 4 is a 6 bit image. The signal to noise ratio for channel k image is defined 


in Decibels as 

(P-P)k 

\ — T“ 


where (p-p), = peak to peak signal value 

tv 

= 127 for k = 1,2,3 and 63 for k = 4. 


In comparing the results of this study with other methods, bit rate vs. 
average mean square error curves and bit rate vs. average S/N ratio curves are 
generated (Figures 1 and 16). This is done by calculating the average mean square error 
for a fixed bit rate as 


2 

e 


% 


4 

Z 

k=l 


2 

®k 


and calculating average S/N ratio for a fixed bit rate as 

4 

s = % ^ \ 

k=l 

These definitions are consistent with the definitions used by TRW^ Inc. in their 
study. The implications shortcomings of these average performance indices will 
be discussed in section 5. 


2. Fast KLT Coding of 255 x 255 Images 

a. Figure (10) shows the 255 x 255 encod.d image of the Channel 1 earth data 
image with 2 bits/pixel average, 1 bits/pixel average, 0.5 bits/pixel 
average, and 0.25 bits/pixel average. Figures (11) through (13) are the 
255 X 255 encoded images for Channel 2, Channel 3, and Channel 4, respectively. 
All these images were photographed on a Dicomed scanner after contrast stre- 
tching each image according to its histogram. 
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FIGURE 10 Fast KLT 255 x 255 encoded image of the Channel 1 Image 
(a) 2 bi*" Sixel (b) 1 bits/pixel 

(c) 0.5 bics/pixc.1 (d) 0.25 bits/pixel 
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FIGURE 12 Fast KLT 255 x 255 endoded image of the Channel 3 Image 
(a) 2 bits/pixel (b) 1 bits/pixel 

(c) 0.5 bits/pixel (d) 0.25 bits/pixel 
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b. Table 2 lists the calculated mean square error and signal (peak to 
peak) to noise (root mean square) ratio of the 235 x 255 encoded image 
for each channel earth data images. 

c. Figure 14 (solid lines) shows the bit rate vs, means square error of 
table (2) . 

d. Figure 15 (solid lines) shows the bit rate vs. signal (peak to peak) 
to noise (root mean square) ratio of Table 2. 

e. Table 2 also lists the average mean square error of each channel earth 

data images. 

f. Figures (1) and (16) show the comparison of these results with the re- 
sults of TRW, InC. effort cited above [20]. 

g. Figure (17) shows the bit allocation pattern of a 255 x 255 image at 1 
bit/pixel average for Channel 1. The value of n^^^ (number of bits for 
the first pixel) is a design parameter, in our experiment, several dif- 
ferent values of n^^^ were tried and the values between 6 to 8 (for p=0.93) 
were found to be most suitable. In order to simplify the quantizer design, 
all the elements in the FKL domain that are assigned equal numbers of bits 
may be assumed to have equal variances. 

h. Figures (18) and (19) show the output histograms of Channels 1 and 3 images 
with 1 bit/pixel average. 

In calculating the average bit rate, we define a parameter (called pseudo bit 
rate constant) PP. For different values of PP the corresponding values of the 
average actual bit rate denoted by ACURAT are printed out. Then the user can pick 
a desired average bit rate, i.e. ACURAT, and from this list [see Figure (20) ] find 
the corresponding value of PP. The value of PP is used in the fast KLT coding al- 
gorithm program. 
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255 X 

255 Encoded 

Images 






J 

Channel 1 

Channel 2 

Channel 3 ’ 

Channel 4 

Average Over All Channels 

Bit Rate 

M.S.E. 

SNR 
in db 

M.S.E. 

SNR 
in db 

M.S.E. 

SNR 
in db 

M.S.E. 

SNR 
in db 

Average 

M.S.E. 

Average 
SNR in db 

2 

2.48 

38.13 

7.39 

33.39 

3.10 

37.19 

0.88 

36.54 

3.46 

36.31 

1 

4.86 

35.20 

14.46 

30.47 

8.11 

32.96 

2.27 

32.43 

7.43 

32.77 

0.5 

10.12 

32.02 

30.05 

27.28 

17.80 

29.57 

5.19 

28.84 

15.79 

29.43 

0.25 

18.02 

29.51 

53.87 

24.76 

29.93 

27.32 

9.13 

26.38 

27.14 

26.99 


TABLE 2 M. S.E.: Mean square error 

SNR in db: Signal*(p-p) to noise (rms) ratio in db 

*: p-p values are 127 for Charnel 1, 2 , 3 

63 for Channel 4 
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FIGURE 17a: Bit allocation pattern for a 

255 X 255 image with 1 bits/pixel 
average. The Figure shows bits 
at every fourth location in the 
upper left quadaint of the trans- 
formed Channel 1 Image. 


FIGURE 17b; Bit allocation pattern in the 
upper left 16 x 16 region of 
the transformed Channel 1 Image. 
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FIGURE 18 HISTOGRAM of the 255 x 255 encoded image 

of the Channel 1 Image with 1 bits/pixel Average 
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FIGURE 19 


HISTOGRAM of the 
of the Channel 3 


255 X 255 encoded image 
Image with 1 biti /pixel Average 
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3 . Faat KLT Coding o£ 15 x 15 Image Blocks 


In. block by block coding, the correlation parameter values as indicated in 
the last two columns of table 1, were used. These are the average horizontal and 
vertical correlation parameter values for a 15 x 15 block. 

a. Figures (21) and (23) show the block by block encoded image at vari- 
ous bit rates of Channels 1, 2, 3 images respectively. 

b. Table 3 lists the calculate-.. .. ..ans square error and signal (root mean square) 
to noise (root mean square) ratio at various bit rates of the four 

block by block encoded images. 

c. Figure (14) (dotted lines) shows the bit rate vs. mean square error of 
Table 3. 

d. Figure (15) (dotted lines) shows the bit rate vs. signal (root mean square) 
to noise (root mean square) ratio of Table 3. 

e. Figures (1) and (16) also compare the average performance of 15 x 15 block 
coding scheme with other methods. 

f. The histograms of these encoded images were very similar to those of 255 x 
255 block size encoding and are therefore not included here. 

g. Figure 24 shows the bit assignments in the transform domain for 15 x 15 
block coding of Channel 1 image. The (15 x 15) block encoded images were 
not contrast enhanced and the corresponding original image is included in 
Figures 21 through 23 for comparison between the original and encoded images. 

4. Other Coding Experiments 

In previous experiments, uniform quantization of x^^ in the transform domain 
was performed. Further, constant correlation parameter values were used in each 
dimension, over all the image blocks. However, the value of the correlation para- 
meters can vary over the different image regions. We examined the effect of these 
variations as well as the effect of different quantization schemes. The various 
experiments are as follows; 
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FIGURE 22 FAST KLT block by block encoded image of the Channel 2 Image 
(a) 2.009 bits/pixel (b) 1.058 bits/pixel 

(c) 0.591 bits/pixel (d) Original Image 
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FIGURE 23 FAST KLT block by block encoded image of the Channel 3 Image 
(a) 2.000 bits/pixel (b) 1.067 bits/pixel 

(c) 0.591 bits/pixel (d) Original Image 
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FIGURE 24 Bit Allocation Pattern for block 
Bit rates 
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ly block (15x15) encoding of Channel 1 image wit; 4 different 



(4-1) Take the average value of horizontal and vertical correlation para- 
meters of the entire 255 x 255 image as a constant correlation parameter 
value to be used in both dimensions. 

(4-2) Use the actual horizontal and vertical calculated correlation para- 
meter values of the entire image. (The image has different values of the 
correlation parameter in each direction) . 

(4-3) Divide the image into 15 x 15 small image blocks. Then calculate the 
vertical and horizontal correlation parameters of each small image block, and 
take the average over all image blocks. The average of these horizontal and 
vertical correlation parameters is taken as a single constant correlation 
parameter value for block by block coding. 

(4-4) Here the average horizontal and vertical correlation parameters of each 
15 X 15 small image block are taken as the horizontal and vertical correlation 
parameter values. 

(4-5) For each 15 x 15 small image block its own horizontal and vertical cor- 
relation parameters are used. We call this 'variable p coding method, or 
'adaptive coding method' since the quantizer for each block will be different. 


The probability density model of the transform domain variables x^^ is 
assumed to be the two sided exponential function 
p(x) = ce , 


ffhere 


xMAX ^ I I 
c = l/( / e “l^'Ux) 

xMIN 


a - 


la 2 

X 


0 2 = variance of x. 

X 


/ 1 

xMAX = maximum value of x^^ 

xMIN = minimum value of x^^ - - xMAX 
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In some of the experiments, a compandor design was used for quantizing the 

A 

variables compandor is implemented as shown in the accompanying 

figure. 



X*. 

ij 


Nonlinear 

function 

The nonlinear f u cion f (x) is given by 


xMAX [ 1 - exp (- ] / [ 1 - exp — ) ] x > 0 


f(x) 


f([x|) , X < 0 


-1 


The output X* is given by f as 


l“riog. z>o 


f ^(z) = X* 


a e 


f-f^(jzl) z<0 

In our experiments the values of xMAX and xMIN were chosen as follows 
xMAX = X mean + 2.5 * SQRT (VAR) 
xMIN ~ X mean - 2.5 * SQRT (VAR) 

where 

X mean = mean value of in the image block 


VAR = variance of 

The value VAR can be obtained either by actual calculation or by the theoreti- 
cal formula for E[x^^] given in Chapter III. 

Results of many experiments that were tried are recorded in Table 5 • 
a. Figure (25) shows the block by block encoded channel 1 image with 
variable p in each block. These images also have not been contrast enhanced. 
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FIGURE 25 fast KLT adaptive block by block encoded image (with variable o in each block)oi 
the Channel 1 Image (a) 2.010 bits/pixcl (b) 1.076 bits/pixel 

(c) 0.616 bits/plxel (d) 0.126 bits/pixel 




b. Table 4 lists the mean square error and signal to noise ratio for adaptive 
coding of Channel 1. Figures (14) and (15) show the performance of adap~ 
tive coder for Channel 1 with respect to 15 x 15 block coding and 255 x 
255 entire image coding. 

c. Table 5 compares the mean square error for 10 different coding experi- 
ments. 

Explanation of Table 5 

c-1 Entire: Treat entire image as a 255 x 255 image block. 

c-2 Block : Divide entire image into 15 x 15 small image blocks then 

perform the coding algorithm. 

c-3 Constant p: Use the average correlation parameter value for all 

image blocks. 

c-4 Variable p : The correlation parameter value varies for different 

image block. 

c-5 Uniform Quantizer: Uniform quantizer in the transform domain is 

used. 

c-6 Compandor: Use compandor in the quantizer design. 

c-7 Actual Variance: Use calculated variance of each image block to 

find xMAX and xMIN. 

c-8 Theoretical Variance; Use theoretical formula to calculate the 
variance of image samples x^^ in the transform domain and then 
find xMAX and xMIN as mentioned before. 
c-9 Bit Rate; Average bit rate for the entire image. 
c-10 M.S. Error: Mean square error 
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Bit Rate 

2,010 

1.076 

O.l 16 

0.126 

Channel 1 

M. S. E. 

1.8534 

5.0819 

8.8548 

i 

22,3115 


SNR in db 

39.3964 

35.0159 

32,6042 

28.5907 


TABLE 4: S/N ratio and M.S.E. vs. Bit Rate for adaptive coding of Channel 1 image. 

M. S.E. ; Mean square error 

SNR in db; Signal (p-p) to noise (rms) ratio in db scale 



Experimental Sequence 


Entire (255 x 255 block) 


ck (15 X 15 block) 


« p_ horizontal 


Constant p 


P vertical 


Variable P 

Uniform Quantizer 


Compandor 


Actual Variance 


: Theoretical Variance 


Bit Rate (bits/pixel 




■ . .1 ... t 

Yes 

Yes 



Yes 

Yes 

Yes 

Yes 


Yes 


Yes 

1 

Yes 


Yes 



1.0044 1.0044 1,0044 1.0044 0.9961 ' 0.9961 0.9961 0.9961 


Error 6.4697 6.0278 5.8556 5.1517 5,5427 5.4579 4.8839 5.2107 

TABLE 5 Mean square error comparison of 10 different coding experiments for the earth data channel 1 image 
’^•'The box marked Yes means applicable, the empty box means not applicable. 



5.0305 ; 14.2570 



















5 . Discussion and Comparsions 


Complexity of the Fast KL Algorithm 

The computational complexity of the fast KL algorithm in terms of number of 

computations is the same as that of cosine, Fourier and hybrid cosine/DPCM coding 

2 

algorithms, i.e. 0(N log^N) computations are required for coding an N x N block 
of the image. The computer memory storage requirements are same as those of any 
two dimensional transform coding algorithm. 

Performance Comparison with Respect Image Block Size 
Figures 1 and 16 show the comparison between 15 x 15 and 255 x 255 image 
block size. See Tables 2 and 4 also. It is clear that w.r.t. both the mean 
square error and S/N ratio 255 x 255 image block size give better compression. 

On the average, the 255 x 255 block size has maximum advantage over 15 x 15 
block size at higher bit rates. For example, at 2 bits/pixel the 255 x 255 block 
size gives 1.5 Db better performance than 15 x 15 block size. This difference 
becomes less significant at lower rates e.g., it is 1 Db at 1 bit/pixel and 
negligible below that rate. This difference is not unexpected because as the 
block size increases, larger redundancy can be removed. 
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Performance Comparison Among Diffe rent Channels 

Figure (14) shows the variation of mean square error of channel image as a 
functlOR of bit rate. Although the absolute m.s.e. values differ siguif ioautly, 
these differenoes are understood by studying the varianoes of the oorresponding 
iuages (see Table 1). M the m.s.e. value of each encoded image (at a fixed bit 
rate) are normalised by the variance of the image in Table 1. then it is seen that 
the normalised m.s.e. are nearly equal for channels 1 and 2 and the same holds 
for channels 3 and 4 also. For example, at 1 blt/rate; for 15 x 15 block Size 
the normalized m.s.e. values for channels 1,2,3 and 4 are .046, 0.051, 0.087, 
0.075 respectively. With respect to m.s.e.. then the four channels may be ranked 
as 1,2,4 and 3. However, on the basis of S/N ratio (see Figure 15 ) the ranking 
is 1,3,4, 2. This is because channels 1,2,3 images are 7 bit/plxel images giving 
a peak to peak value of 127 and channel 4 images are 6 blt/pixel giving a peak 
to peak value of 63. If the S/N ratio were defined in terms of the ratio of the 
variance of the original signal to the m.s.e., then the two criteria will give 
a consistent ranking of channels. 
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Comparisons With Other Coding Methods 

Perhaps the most interesting (and surprising) results of this study are the 
comparisons of the fast KL coding algorithm with other coding experiments per- 
formed by TRW, Inc. (see [20]) on the same data. Figure 1 shows the performance 
of the fast KL transform coding algorithms vs. other methods. It is clear that 
the 2 dimensional fast KL algorithm (255x255) as well as (15x15) dominate the 
other methods including the 3 dimensional methods viz KL (temporal domain) - 
Cosine/DPCM and KL-Hadamard/DPCM. However, these curves have to be observed with 
certain caution. First, the mean square error is the average mean square error 
over the four channel images. A more meaningful comparison could be to compare 
the average of the normalized m.s.e. of each channel, where the normalized m.s.e. 
is the ratio of m.s.e. of a channel image and the variance of the original image 
signal. Also, the various methods should be compared on each individual chan- 
nel as well. Unfortunately, the m.s.e. (due to other methods) per channel is 
not available in [SO] and hence these comparisons cannot be reported here. 

Figure 16 shows the comparisons with respect to average S/N ratio. The 
255x255 fast KL coding scheme performs better than the other schemes. The 15 x 
15 block fast KL coding performs better than 2-dimensional cosine as well as Cosine/ 
DPCM at 1 bit pixel by more than 1 Db. At higher bit rate, e.g. 2 bits/pixel, the 
difference between 15x15 fast KL and 2-dimensional cosine is as much as 2,5 Db. 

This difference is quite surprising. The cosine/DPCM and 15x15 fast KL have nearly 
equal performance at 2 bits /pixel. 

Ihe difference in the performance of fast KL vs. Cosine may partly be attri- 
buted to the different correlation parameter values used. For 15x15 block coding 
we have used lower values (see Table 1, p .6 to .7) whereas in the TRW study 
p was taken ^ .9. This difference will cause the bit assignment in the transform 

domain to be different. 


65 


L 


Another interesting comparison is noted v?hen Figures 1 j.nd 16 are compared. 

While 15x15 fast KL gives lower average mean square error (at 2 bits/pixel) than 

the 3-dimensional schemes, it does not give a higher average S/N ratio, This is 

because averaging the m, s.e, over all channels does not average the S/N ratio 

over these channels. „ .r,r,r 

REPRODUCIBILITY Or’ ' 

ORIGINAL PAGE IS POi)U 


Adaptive Coding 

In the adaptive coding method we changed the bit assignment from block to 
block according to the variations in the correlation parameters, keeping the 
average bit rate of each block fixed. The improvement of this adaptive coding 
scheme is only marginal over the non adaptive fast KL. However, a better per- 
formance can be expected if the average hit rate per block is varied according 
to its variance and the bit assignment within each block is based according to 
the correlation parameters of that block. 


6 , Conclusion 

1. To the extent of the usefulness of the average mean square error and 
average S/N ratio criteria, the fast KL transform coding algorithm 
gives superior performance as compared with other 2-dimensional coding 
methods for multispectral images. 

■2. More detailed comparisons are required to establish the usefulness of 
this algorithm from a user’s point of view. These comparisons should 
include other criteria such as (i) the spatial distribution of errors 
in the reconstructed signal (ii) Classification accuracy based on 
ground truth. Also a larger set of images should be used to vali- 
date the conclusions. 

3, The fact that correlation parameter values change sign! '’leant ly as the 
image block size is changed suggests non-stationarities in the images 
and fast adaptive algorithms could be useful. Since, the fast KL 
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transform basis vectors are invariant to changes in image statistics, and 
the transform domain variances are known in closed form, application of 
this method for adaptive coding is suggested. 

4. Since the Cosine/DPCM algorithm has been shown to perform significantly 
better than 2-dimensional cosine [20] at high bit rates, an interesting 
question is, "How does fast KL/DPCM perform"? 


67 


REFERENCES 


1 . 


2 . 


3. 


4 - 


5. 


6 . 


7. 


8 . 


9. 


10 . 


11 . 


12 . 


13. 


14. 


A. K. Jain, "A Fast Karhunen-Loeve Transform for Finite Discrete Images,” 
Proc. National Electronics Conference , Chica V), Illinois, October 1974. 

H. C. Andrews, Computer Techniques in Image Processing , Academic Press, 
New York, 1970. 

W. H. Chen, "Slant Transform Image Coding," Ph.D. Thesis, University of 
Southern California, May 1973. 


W. K. Pratt, "A Comparison of Digital Image Transforms, The Mervin J- 
Kelly Communication Conference, Univ. of Missouri, Rolla, October 1970. 

A. Habibi and P. A. Wintz, "Image Coding by Linear Transformation and 
Block Quantization," IEEE Transactions on Communicat ion Technology.* 

Vol. Com-19, No. 1, February 1971. 

N. Ahmed, T, Natarajan, K. R. Rao, 'Discrete Cosine Transform', IEEE 
Trans. Comp ., Vol. C-23, No. 1, January 1974, pp. 90-93. 

P. A. Wintz, "Transform Picture Coding," Proc. IEEE , Vol. 60, No. 7, 
pp. 809-829, July 1972. 

W. K. Pratt, "Karhunen-Loeve Transform Coding of Images," P roc. 1970 
IEEE Int. Symp. Inform. Theory , 1970. 

H. P. Kramer and M. V. Mathev/s, "A Linear Coding for Transmitting a Set 
of Correlated Signals," IRE Trans. Informati on Theory, Vol. IT-2, pp. 
41-46, September 1956. 


H. Hotelling, "Analysis of a Complex of Statistical 
Principal Components," J. Educ. Psychology, Vol. 24, 


Variables Into 

pp. 417-441, 498- 


520, 1933. 


W. D. Ray and R. M. Driver, "Further Decomposition of the Karhunen Loeve 

Series Representation of a Stationary Random Process, 

mation Theory , Vol. IT-16, No. 6, pp. 663-668, November 1970. 


V. R. Algazi and D. J. Sakrison, "On the Optimality of the Karhunen-Loeve 
Expansion," IEEE Trans. Information Theory (corresp.) pp. 319-321. 


S. Watanabe, "Karhunen-Loeve Expansion and Factor Analysis, Theoretical 
Remarks and Applications," Trans. Fourt h Pra gue Conf. Inform. Thgp rx, 
Sta tist. Decision Functions and Random Processes , Prague, iybO, pp. 

635-660. 


A. K. Jain, "Image Coding via a Nearest Neighbors Image Model," IEEE Tran s 
Communications, Vol. COM-23, Nc. 3, March 1975, pp. 318-331. 


68 


15. A. K. Jain and E. Angel, "Image Restoration, Modelling and Reducation of 
Dimensionality," IEEE Trans . on Computers , Vol. C-23, No. 5, May 1974. 

16. I. J. Good, "The Interaction Algorithm and Practical Fourier Analysis," 
J. Royal Statistical Soc . , (London) B20, 361 (1958) . 


17. E. 0. Brigham; The Fast Fourier Transform , Prentice Hall, New Jersey, 1974. 

18. A. K. Jain, "A Fast Karhunen-Loeve Transform for Recursive Filtering of 
Images Corrupted by White and Colored Noise," (to appear). 

19. A.K. Jain, "Image Restoration by Operator Factorization of Point Spread 
Function and Image Covariance Function," Dept, of Electrical Engineering 
SUNYAB Report 75-001, Oct. 1975. 

20. A.Habibi, Study of On-Board Compression of Earth Resources Data, Final 
Report, Contract No. NAS2-8394, TRW No. 26566, TRW Systems Group, Redondo 
Beach, Calif. 90278, Sept. 1975. 

21. J.W. Woods, "Tt^o Dimensional Discrete Markovian Fields," IEEE Trans. 
Information Theory. Vol. IT-18, No. 2, March 1972, pp. 232-240. 


69 


APPENDIX I 


CORRELATION PROPERTIES OF THE MINIMUM VARIANCE REPRESENTATION 


The equations of representation of the first order stationary sequence {xJ 
with zero mean and autocorrelation 

E[x.x. ] = p ^ ^ ^ ^ (1) 

1 2 


are given by 


X. = ct(x,,- + X. .) + V. i = 1,...,N 

i 1+1 1-1 1 


^0 " ^^1 ''' ^ N +1 ^ N +1 

where a = p/(l+p^). Multiplying (2) by taking expectations and using 

(1) we get 


( 2 ) 

(3) 


pl^l = aCpi’^ + E[v^x^^^] 

2 

which for jk| - 1 and with a = p/(l+P ), gives 


E[v.x. .,] = 0, kj^O l^i^N 
1 i+k 


(A) 


(5) 


2. A „2. 


Similarly for k = 0, E[x^v^] = (l-p'^) / (1+p'^) - 8 ; 1 A i < n (6) 

Multiplication of (2) by taking expectations and use of (5) and (6) yields 

^^Vi+k^ " ^k,o 




(7) 


where 6 is the Kronecker delta function. Similar procedure when applied to 

ij 

equations (3) gives 


and 


E(V(,x^l = (i-pb«^ „ ; E[Vk' = 




( 8 ) 

k,N 

( 9 ) 


I-l 



APPENDIX II 


A TWO DIMENSIONAL IMAGE REPRESENTATION 


Let {u. . i, j = 0,1, . . . ,N,N+1} be a two dimensional stationary sequence 
1 j > 


with zero mean and a separable autocorrelation function given by 


( 1 ) 


Let u be the best linear, mean square, estimate of u.. obtained from all 
13 ij 

u not including the point u... This is obtained by first writing 
mn ° ^ Lj 


“ij ^ ^ \)i^“i+k,j+£^“i-k, j-t’*'“i-k,j-£^“i-k, j+J.^ 

k+MO 


( 2 ) 


and minimizing E[u^^-u^^) ] over the coefficients Differentiation of this 

expression with respect to a , setting it equal to zero and after some alge- 
braic manipulations we get for 1 - {i,j) - N 


]k| , jk+k’l ^ |k-k'L, |t+si’L 

p' 'p‘ - I j; a^,j^, (p' ' +p' '+p' ') 


k’ I' 


(3) 


for all k + )i 0, k' + J.’ 7 ^ 0. From (3) it can be proven a^^ = p/(l+p ) = 

n 2 2 

a ; a = a^ and a = 0 for k +11 - 4, so that the two dimensional rep- 

Xu XX UX 

resentation equation defined as u.. = u.. + v.., becomes 

ij ij 

u. . = a(u ,+u. ..,+u. . +u. ) 

ij i+l»J i.J+l i-l.J 


- a (u.,, . 1 +u.,, ., 1 +u . ,) + V.. 

1+1, 3-1 1+1, 3+1 1-1, 3+1 i-l,J-l ij 


(A) 


where a = p/d+p*^). Eqn. (A) can be rearranged to give 


u^ . - ot(u. , , .+u. .) = V. . 

i3 1+1, J 1-1, J ij 


V.. - a(v. ..,+v. , ,) = V.. for 1 - (i,j) - N 
ij 1,1+1 3,3-1 11 


(5) 

( 6 ) 


Following Appendix I, the statistical properties of v.^ and are obtained, 
for 1 - (i,j) - N, as 
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t 


E[v..l = 0 = E[v..] 

“‘l-l.k’ 

> = 8'‘<“ik-“‘i+i.k-“'‘i-i.k’'*3r“'3+i.*'“‘3-i.i‘' 


where 


(7) 

( 8 ) 
(9) 


8^ = (^) . 

1+p 

At boundary points (i=0,N+l; j=0,N+l) the minimum variance representation equa- 
tions are somewhat different. These equations are not given since they are 
not required for the derivation of the fast KL transform result. Similar 
results are obtained when ^ above. 


RKPT?OTlUCi.nTF.ITY Or ■. 
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APPENDIX III 


USAGE OF COMPUTER PROGRAMS 

Two computer programs written in FORTRAN IV have been developed. These are: 

a. image ANALYSIS PROGRAM 
h. FAST KL TRANS ORM CODING PROGRAM 

These programs have been implemented on the NASA/MSFC IBM 360/65 and NASA/Ames 
CDC 7600 computers. This chapter explains the usage of these programs. 

1. IMAGE ANALYSIS PROGRAM 

The MAIN routine of this part of computer program does the following 

jobs: 

a. Print out the histogram and statistical parameters of the image 

b. Calculate the horizontal and vertical correlation parameters over 
all image blocks 

c. List a pseudo bit rate constant versus actual bit rate average table 
The subroutines used in this program are: 

a. KAR - This subroutine finds the maximum, minimum, and mean for each image 

block, creates zero mean image block and calculates Pq-|^ 

and bit allocation for each pixel. 

b. HISTO-Prints out the maximum, minimum, mean, standard deviation and 

the histogram of the 255 x 255 image. 

c. CORRE- Calculates the horizontal and vertical correlation parameters 

of 255 X 255 image. 

The definition of all the parameters used in these subroutines have 
been clarified in the computer listings. 

2. PROGRAM USER'S GUIDE OF IMAGE ANALYSIS PROGRAM 

A 256 X 256 image is loaded on a standard magnetic tape. In order to 
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process via the Fast ICL Transform technique, the actual dimension of the image 
used is 255 x 255, because if N is the image size N+1 must be a power of 2 to 
implement the sine transfom. 

2. 1 PROGRAM INPUT 

The JCL cards of the image are shown here as a reference for the user. 
//GO.FTOlFOOl DD DSN=INPTl, UNIT=TAPE 9, VOLUME= (PRIVATE, ,SER=A0001) , 

// DCB=(DEN=2, LRECL=512, BLKSIZE=516, RECFM=VS) , LABEL=(1, BLP,, IN), 

// DISP=(OLD, DELETE) 

This means the image is loaded in the following format 

a. 9- track magnetic tape 

b. 800 BPI (bits per ir. h) 

c. Integer* 2 (2 logical ds to represent one integer pixel data) 

There is one input data card and it is; 

$INPUT IFRAME=1, N=15, Nll=7, ILINE=17, $END 

2.2 PROGRAM OUTPUT 

Output listing of this program includes 

a. Histogram of 255 x 255 entire image 

b. Horizontial and vertical correlation parameters of 255 x 255 entire 
image 

c. Print out a table of pp(I) and ACURAT(I) , which allows the user to 
pick any desired bit rate of the image to be used later in the Fast 
KL Transform coding program. 
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3. FAST KL TRANSFORM CODING PROGRAM 


This computer program does the following jobs. 

a. Creates differential image of a 15 x 15 image block, 

b. Applies Fast KL Transform to 

c. Calculates bit assignments to different elements in the transform domain. 
Performs auantization. 

e* Applies inverse Fast KL Transform. 

Stores final result as a 255 x 255 image on magnetic tape. 

S* Prints out the Histogram and the results of analysis of the encoded 

image. 

The subroutines are: 

a. KAR -Creates a zero mean image and computes bit allocation pattern 

of a 15 X 15 image block. This subroutine is linked to sub- 
routines CODE, qUANT, RECON. 

b. CODE-Creates a differential image of a 15 x 15 image block, calculate 

the mean, variance, standard deviation of the transform domain 

15 x 15 image block. This subroutine also calls subroutine XFORM. 
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REPRODUCIBILITY OF THE 
ORIGINAI, PAGE IS POOli 

The definitions of IFEAME, N and Nil are given in program listing, parameters 
PAIOj PAOI, and PCON are obtained form the Image Analysis Program. PCON har 
a same meaning as PP, the pseudo bit rate. In our example PCON = 0.125, 
was picked from PP(96) , and the corresponding value of ACURAT (96) = 0.99607837 
which almost equals the one bit/pixel average rate. 

In this program, two magnetic tapes are used, one for the original 255 x 
255 image, the other for storing the final 255 x 255 encoded image. The JCL 
format cards are shoiim here for the user’s reference. 

Tape 3 (original 255 x 255 image) 

//GO.FT03F001 DD DSN=INPT1, UNIT=TAPE 9, V0LUME= (PRIVATE, ,SER=A0001) , 

// DCB=(DEN=2, LRECL=512, BLKSIZE=516, RECFM=VS) , LABEL=(1, BLP,,IN), 

Tape 4 (final 255 x 255 encoded image) 

//GO.FT04F001 DD DSN=0UPT1, UNIT=TAPE 9, VOLUME=(PRIVATE, ,SER=SAVE) , 

// DCB=(DEN=2, LRECL=1024, BLKSIZE=1028, RECFM=VS) , LABEL=(1, BLP,,0UT), 

// DISP==(NEW, DELETE) 

c. QUANT -Per forms quantization on the transform domain 15 x 15 block 

image samples 

d. XFORM- Initializes and rearranges the data to call subroutine HARM 

and implements the Fast KL Transform on a 15 x 15 image block. 

e. RECON -Reconstructs the 15 x 15 encoded image block from transform 

domain quantized 15 x 15 image block. 

f. HARM--An FFT subroutine which is available in the IBM scientific 

subroutine package. 

g. HISTO-Calculates the statistics and mean square error of 255 x 255 

final encoded image and prints out the output histogram of the 


same image. 
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4. PROGRAM USER*S GUIDE OF FAST KL TRANSFORM CODING PROGRAM 

4.1 PROGRAM INPUT 

First select the actual or desired bit rate for encoding the image (say 
1 blt/pixel). Then the Input data card for 1st channel image, 15 x 15 block 
coding, and horizontal and vertical correlation values of 0.^1, 0.598 respoc- 
tively is; 

$INPUT IFRAME»1, N»15, Nll-7, PA10»0. 71059489, PAOl-O. 59758818, PCON-0.125, $END 

4.2 PROGRAM OUTPUT 

Output listing of this progrcm is simple and includes 

a. The actual bit rate of 255 x 255 encoded image. 

b. Output histogram of 255 x 255 encoded image, various statistical 
parameter, the mean square error in encoding, and the signal to 
noise ratio of the encoded image. 

The listing of subroutine HARM is not given in this report, since this is a 
standard program available in the IBM scientific subroutine package. 
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APPENDIX Vt 


IMAGE ANALYSIS PROGRAM 


C 

C 

c 

r 

c 

r 

C 

C 

G 

C 

C 

C 

C 

c 

c 


THIS p^OGkAM uDtS FOLlOhIIIG JOSS* 
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?. CALLULATE horizontal AMD VERTICAL GGRKELAT lOK FAkAMETERS OF 
ENTIRE IHA G£ 

3, FTM.i iOiTrtiin. FSEUuO UlT FaTF CONSTANT FcR DIFFERF-NT BIT 
tvATt. ASSiGUtIfNT 


lU : 
Fl'J ! 
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HAi: : 

FmTI » 

ACURmT! 
IJ ! 
s 

« M R s 
FP ! 
vMwX t 
XHIN : 
XfiuAN : 
XVAR 
rA NCl 
1 G TmL 
Su 
F 1 ‘J 
F «1 


hORIZONTAl. INPUT uATA STRING 

HOrIZC.UAL GORRcLATICN FARAiTETER OF AN IMAGE BLOCK 
Vi='kTiL.Aw t;0KR]Li.ATIuM PARAHETCK GF AN IMmOc HtOuK 
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NAMELiSTx INFUT/IFR AME ,N*Nll, ILINE 
NAME^I jT/OuTFUT/FAlu* PATl 
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ALIN--Il 1N‘P 
At^sN 

GONST-17 .* AtiNE 
DO Eo I=^l»16 
DO 2- lo 

u ti» n-c . 


EEPRODUCIBILITY OF 'HIE 
ORIGINATE PA<1R IS POOR 


. W W WM»» WWIp 


i MWa i Mllli 


?3 rOMTiMUE 

DO 21 
ACK(II=C« 
21 rONTlMUt, 
t-oii.) s ?.5 


DO 2 2 I=2,4«'? 

ruNTlHUt 

on 1 . IC= 1 ,ILIN£ 

CO 11 J=1 ,15 

PtflCX*) UUIK) .<=^i,25b» 

00 12 1=1 rCbS 

10(1, Jl=i JSI* 

12 CONTr’JUE 
11 C;^NT1 4 0£ 

0^ iO 1=1,17 
iS=(I-U-l5 + l 

DO IV J=*,15 

ro lu ii=rs,isp 
II=x\-IS+l 
U4U, j» = It2UA,J) 

14 CuNflNUc 
CaL*. tfflii 

13 rONTlUUt 
ID CO.UlJUt 

Ffllu = '’4iC/C0Mr 
>AC*=!>ACi/rcMsr 
UUfi=H‘J>AN*CCNST 
TOTAL=XMEAN 
Xf1tA'4=Xi'F AN' /l)L» 

X tf An= n 1/ Mk / COM - X^it ft M“ X rt E C M 

00 3 - Jrt = l , Si. C 
. ACU^Mr {JA > = mCR(JA ) 70UM 
5D CONTINUE 
RE'/UNO 1 
C.AU HISTO 
RPt^l'n 1 
UAU COPKC 
WRITE (5,^11) 
irtPm.'otOUTPUT) 

00 21 1=1, 4DG 

.CUTE (orllg) I, Udi , I, ACURAT(I) 
25 r,CNlI'.'Ut 
STOP 
‘•'NO 


IV- 2 


BffiKODD-caBttnT OF Tm 

OmnSAL PAGE is P0,T 


KAP> 

CuMil JrvANALYS/XKAX.XMlMt XMEANf XVAR » TOTAL 

OU.l:ini-</okOUP/UTiD, ibi iPPIAOCl »NtNll«MuMHDAIl5 il5) 
rUMMOl/THI»«G/PAlCi‘*A01tACBi(4CC) 

SU'n=3. 

SUM? = ^, 

on i: j=i.N 
ou i: 1=1, N 
suiU-.Tam^uu , j) 

SiJMc; = L»iJM? fU ti , J)*U II , .1) 

iFfJd.J) .GT. xr.MX) XMAX = U(I,J^ 

*r(’J(i,J) .LT, XiliN> XiilN=U (I , J> 

10 CoHTiNUF 

A»(£A J = '5UM1/ (N*N) 
xiiE;*,j=xM^ftrHUurn 
X \l >*R= A \l tik. SUi'i I 
SUM=0. 

00 1- J=1,N 
00 11 1 = 1 ,N 
0 n , J) =U (T , J) -AMLftN 
SU.1=j’JM + U U ,J) 

11 OCNTiNUt: 

VA*^=3UM/ (M»N> 


21 


P1Q = ^ . 

PjUJ. 

OL 21 J=1,N 
DO 1 = 1, N 

P13=310*U(1,J)*U(I«-1,J) 
1="''? i + lMI , Jl *0 II , J* 1) 
CC'iriNUt 
P^ii=^lG/ I 
PCl= 1/ t JAR^N*M» 

PAlC = i- AiC *F1 C 
PA0l=rAOl*-FOl 
AuPlsPlO/d.^-Flw^Pii) 
AlD 2 = P jl/ iI.+Pjl*POl) 


ANSJ = 'J*N 

Am=o. 141592F/ rN*l» 


suM=: . 

DO s’ 1=1, N 
DO 5 J J= 1 , M 

DoM- 11 .-d .♦*.L>l*rOS( I‘AA» »»( 1.-2. 
AwAMOi a, J»=DUM 

i,iJM= lj,i*ALQDlC<0OM/AL AMOA II, 1) J 


*ALP2»C0SU*AAI ) 


DO COMTIMOE 


IV-3 


54 

*53 


CJENs jUM 

an jA=i,<tca 

f-cON^FPCJAl 

nuMsANir.*Ull-PCuu)/OtM 

N£ai=" 

00 J-1,M 

OJ 5<* 1=1, N 
AF.rUL4*4DA tl,j| 

W&0U t = Al4.-iiU{!‘H^Llf.lijl AR/AtAMDAII,!) ) fru ,5 
If (,'iAuUM ,Gr. NA0Ui1=Nll 

iF<NADU!-i ,LT. 0) nA0UV=C 

uSUtti-NAObM^-fiSUM 

CLNlI'JUf. 

MCRt Jh > = AbR ( JA » *u£U4 

CONTI'^UL 

E ^'0 


IV-4 


r WMyViYniitM 




wSSflJ!:S act 16»1 . IEHP.65. . S JMBOU UCC I . I6RAPH ( 10 r.fl 

Co'1hJS/ANALYS/Xf^A<*K*1iN»)(r.tftKi)fi<»*R.TOTAL 

LoMrtOM/i*P«Ct/ lU 12561 

INT'G=R»2 lu IH^/.'SlANK/IH /* 


l"jA‘^rr^/6i-u/.TtllP/65*3./.J6SAPH/lJ*!/.SIAP/lHV.BLANK/lH /, 

•uA5H/.H-/,FU»S/lH>/ 


Fijk^vT t/y///) 
FuRMST C Jrt'?> 


l-ORMAT CD'tlSafll) 

F C 5 T ( 3 X f 12 1 i£ A f t A 3» 5 1 3 X * Ib 1 1 X j 1 M '* # 1 0 0 A 1 1 

RA''IGt:= XKftX-XrtiM 
S a=SO<T tXi/Aft) 

W R I T_ (6flJ*!)' 

rt P IT- ( 6 f i. 3 ■» > >* tiA X f XM I u t R A N6E , Ti- < /*»- * XHE AM t '^C> 
XlMt= t xwfix-XttIf4i /bA. 

TEMF <1 )sX4lN 

no X 1=1 f’y*^ 

RuAUtU i rut J > t J*ltt 5 p) 


DC 1 * J=l» 2 ‘?t 
XCJJ=IU{J> 

cOMlIlUT 

DO *♦ .1 = 1*255 

rn x=i* 5 A 

TtMP =Tt>'F IK > f Xisn 

it-(XU)-TtMMK*in b»6*5 

CCNIIUUF 

in CK* =irT|K) ♦! 

rOMiriUF 

CGMTlNUr 

IuT“AX=0 

DC 7 1=1 

iF(lGTtn.GT.ICrMAX) ICTMAX = ICrn» 

uONTlMUt 

Mu-I-^TMAX/AOt 

WRIT-{6,1''1) 

WKlTZ{otia“) ^'0 
to i.’ I=itlC 

I DriA'>rt tn =1*MU-*1C 

rCNTIMUE 
rtRITtt6*l31) 
wRif tbf 1C5) IGF A"H 
DO 1‘* t=ltl32 

tt«Rtli=Of 5H 

CtrtTluUE 


ttEPEGDUCmiOTT p: 

ORTGTNy^L ^ 



00 l-> I = ,iLtl3Ltir 
OARtUspLUj, 

Ifi rONTiKUt 

WRIT";6»1S6J FiAR 
DO <i 1 = 1, 6«» 

■ IFtICTCII .Gt.NOI 50 TO 9 
no 1'' K=1,10C 
sr'iGuwtK* =n*.AnK 
IQ CONTI fJUC 

GO TO 11 

9 J=luT<n/NO 

^ 1 =J ♦■1 
CO li K=1,J 
S (K) =STAP 

12 rouTiNut 

on ij K=Jt,lGt 
S YMii-L (M =eLANK 
15 GCiNTI‘'*Ufc 

WRIT : J 6» 177 I I » T tMr*( aI , Ic T ( i ) , (SYMtlOL (K) , K=l. 13 QJ 

a roWTlMUc *♦ 

NRITuT6,1C6) bap 

(6,175) lOPAPM 

R Fluor. 

tMO 



Suj^iuTir:^ rOnKfi 
□ IMlNSION a <2a6» ,n |£56» 

COM!iaH/*.NALYS/AMAX,XMTN,Xf1tAN,X^A^,TCTflU 
UoMMUN/bPAUE/ 1 UI 2a6) 

Ir4Tc.r.^(<*2 lu 
DATA 3/?*:b*C./ 

120 f uRMAT (/////) 

Itl FORm:;T tlX»»HOKXXowTiiL DOkRtL>»IION PAkAMFTEr OF £MIRt IMAGE -r, 

OF ENTIRE IMAGE =t,F12,0) 

*• 1 j - «i » 

uO n Isi,255 
RFAu(l) 2U 
00 5; Ksi»25t< 

M(KlslU(r» . 

5»? CJNTIMUF 

DO 2: .J=l,254 

Ki *3= (i. I Jj -XMJAW) ♦ (A{J*-*.|»XNEAN**Rn 
2Q ’"O.NTIMUL 

IF (I ,G3 . i> GO rn (* 

Go n I 

4 on --I 

R j 1-{ J { J) -Xht A m* - lACJJ •XMoANi ♦40 I 
ft C^NTlMUr 

7 GO 53 J-i , 2?o 
BJ J» = A ( Ji 
31 r'CNTiMUt 
oONTINUE 

K10 = P*0/ CX»/AR*2i>5.*25r>. > 

F'‘l-^:i/ (XVAR*2‘i5.*2b'1,l 
WRj.TE;(t>,123) 

»JPlT£(b,l 21) Kiw,R0i 
KE TJ^?M 



APPEOT)IX V 

FAST KARHUNEN LOEVE TRANSFORM COMPUTATION TWISTING 

t 

Inia r^KOGKMM DOES FOu*.OMING JOBS! 

1, b-T niFK£w,ENTIAL IMAoE 
£. F-;-<FOkM FASI K. l. XFORM 

3. OMLCOLAlt dir RATE ASilGUMt^Nr 
< USF rUNSTAtiTiRHGi Iw ORuc-R TO GET SAME QIT 
ASSiGuMtNI RATTEf^N FOn EACH IMAGE HLOCK ) 

4. Ar>PL» UNIFCKM QUANTISER rtCHNIOUE 

5. REF FORM FAST K. . Xf ORM 

f, STORE final RFSUlT iNAGt Urt J A»»F 

7. P'^INI OUT COTfUT mNmLTSIS hUC rilSTOGRAR 

iU : HOkIZoNTmL INPUI OAfA STRING 
rCOK » FStUDO oil RuTt. lONSTAMT 

PhI: i i.Vt.R«Gt Or MQitlZOMA*. wORRE«ATIt)N PAivAl'ETtR OVER A^w 

liAGE BLOCKS 

r*A''L I mvlaAge cf vertilml correlation parameter over AtL 

IFhGc 3..orKS 

AIo t 256 JY I.IhGE BLOCK 

N*1 » ASSUMLut liiT RATE FUR T«E MOST SIGNlFioANl PIXtt IN XFORM 

uCMAINr This VAwUE is OMwLUtATEu FPC.1 THtOREflCAu FORMUwA 
mCURAT I AoTUAL AvERAGt BIT RAU- 
AriE«H t I'.fAN LF CRiGiNAL IMAGE CLOCK 
vmR J VARlAVoE UF OrilGIMl*^ IMAGE duOCl^ 

NlEV * 15 tv lb "IT HATTEkN mRRiY 

XMc.h'I t MEAN Or KrORM OuFiAIM IMAOi d^-OCK 

XVAR I VArIANlP JF XruRM OOMAIN IiIm&F oLOCK 

XuEVi t OTAMuARL L'EVIATIOm JF XFJPM UOKAlN IMAGE B^OCK 

fMAX t MAXlHUrl ni^IPF'IKG LExEl uF XrOPK iTOMAlN IMAGE CLOCK 

f Mif. t NiNA‘»UM TRIPPING QP XFoRM CJFAIN IMAGE 0*,0CK 

X.1AX t MttXllUM v/ALUE or FINAL «cSULT If-AGt 

xalM J MINIMUM XA^Ul. of FINA*. kESUvT IMAGE 

PmNl.- t M^xI^.LM tfALUr MINUS MlNlMuM VALUE 

T«TAL » SUM oF TCTA^ VAwUtS 

HfAN OF FIMAl. RESULT IMAGE 
E'U * STANuA^U liEvTATIUN OF FINAL RESULT IMAGE 
F^Rl : MuAN SnUAPt ERROR 

LRP.2 * koOT F.cmN square EkkOK 

tRRuHT MlAM CCOARE EkPOR IN uT SCALE 
LkRoi2» w.wOT SOUAfvb EkRCk IM Qd SGCvE 

£Ni t signal to NOiSt. RATIO t RANGF VS ERR2 > 

Sw? » SiGNAw TU NOI^E kATlO ( SO VS ERR2 I 

SHOP; « SNl IN Ui SCALE 

SNOB’ t SN2 IN Db SCft.-E 

CIMlNCIuN «IB4256j15I 
COMMuN/AMALYS/lU4t5b) »Xl25e) 

OUMM TU/uKuUP/UIl?* 171 ♦PvUN»Nlltf,tPA11»PA>,if ALFltALP2f 
•NlEVI 15,15) ,AL*MLAI15,15I , VAF.Vtlb.i^* »XH4l5t 15> I 
•XmoJ^i., 15) .LiSlAF il3,l*^» iUSTAkH4lb,15) 
r JMMU*1/R|.1 £/NA*A,iaQW 

NAM“LIST/tNHUT/lFKrtME ,N,Nlt ,FAlg,RA5l, ROOF) 

MAMEL rSI/CLilPUf/ACURAT 

INTcG“k*S IU 
REAL*4 Aid 
REAL* 4 X 










ASU=‘»55*2ijR 

KftCif?*'’ 

ifiiiw=': 

L* T ‘4 j J • ^ 1 1 7 
Du 4. 1-1,17 


U(I,JI=U. 

40 rONTI'JUt 

Rf:5u»‘",lfl°UT) 

WkIT C C q , I^4H UT ) 

on » 10=1,17 

on 11 .J=1,15 

Pl 4D<?> ( lull), 1=1, 25 u) 

on Z I=t,25F 

Aioir, j) = iuii ) 

2 ruNTIMUE 

11 rONJiSUE 

00 ea M*i=i,i7 
IS=4.*-15* (fUI-1 ) 

ISP=I^+lt, 
or, i: j=i,ie 
jj=j*i 

00 A=14,ISP 

II=I-IS>2 

ijdi, Jui=Aioa,j) 

20 CCiilMUE 
rouTir-iUd 
IRJh=IB'JW*1 
TALw KAR 

no F<: j=i,i‘i 

DO o'' i = IS,Ii;E 
11=1-1541 

AIdU , J)=USTAkHUI,J) 
fa 5 CuNTT,jUF 
•*? uCNllNUE 

00 6i Jrl,15 
A lBli.vb, J) = C» 

lAIf’ll, J) ,I=l,?5b) 

»v» roNiif'Ut 
5 TuNTINUl 

on 1=1, 25b 
AT01T,1I=C. 

5« CONT IflUP 

WRITl (4) lAltill,!) ,I = x,?5fal 

l,Ch’=') Aufc 

ACUR..'VT = AbK/U6.*i5.) 

r^l-ITE ;b,CUiruT) 

RdWl‘10 3 
R«^rtriO 4 

call histo 
ST DR 



i,u3^0UriNi KAR 

CnMM0U^r.P3lJP/U<17» 171 ,PC0M,M*1|N*PA1C, {^Ai It ALU *ALP2» 
•PL Erf f I5ii0) tALAKUA t VAR tV{i'>*15) ,XH (15tl5) i 

•XHiUStlS) ,U;>TAK(;5ii5l fUSlA RHUat 15) 

CO.IMJN/JUT/Aitt AN 
COMMJ.t/RATE/NALP ti"t»>’) 

NPl=*t + l 
N32sv(t2 
SUM- 3, 

CO i: J=2,NH 
DC I=2fNPl 
SU?1= >'JM*u tliJ) 

10 cuNTiMur 

A‘-£A.)=3UM/ LN*I!) 
sijM=: . 

Cn i.1 J = ?.NP1 
OJ 11 I-?tUPl 
Uilt j) -U4lt J) -AMEAM 


SUrt=SUM«-U(ltJ)*U(I , J» 

11 rcMl INUE 

i/Ar<s:jUii/ (N*N) 

ML»i=P AlC/ 11. ♦■PMit ^PAin) 

ALp2 = '^Ai.1/U. ♦PaC1*PAu1 ) 

AN3U=N*N , 

nA-i. (NH) 

if .oT. ii Ou TO 3U 

SUM= ; , 

era’ I-lfM 

cn sc J=l.N 

DUM= ( 1 AU'i-COSa*AA| ) * U.-2.* ALP2*Cn3ij'^AA>) 

ALAUJA li,J)=:DUM 

suM=su:i*AwOGi(;au'i/A».AMaA(i,in 
*^0 CuNl lMUf 
DENr 3UH 
All=U-.l 

CUH- A^ sn* ( Aii-PCi-NMOEN 
DO 54 jsltN 
DO 5^. l = ltN 
ARsALAMDA (ItJI 

NAuU.l=All-bUT*AwOGlf; I AR/M..MMCAIi,l » > *3 .5 

IFI'UaUP .oT. illi) nA!"jU-I=U* 1 
IFiNHLUM .IT. u) 4.AOUM-3 
NAC-^ = NAOUM*NALk 
NLiW 1 1 1 Jl =NACUM 
54 CuWTlMUE 

JO call conr 

GALL OUAMT 
lhLL Kirnu 
K^TUOM 
fJO 






J! 

I 



i*ua«,UJTlMt. CCO£ 

CUMM0N/GRulJP/Ull7, i7> tPGUW,Ml,U»PAlCiPft Cl.A*.Pi»A>.P2t 
nit.£V( iSfl'i) tALAKuA(15,lbl t tfAP,tf<l5,l5l ,KM.i 15» l5Jt 
“XH3Vl» tl5l tUSTAF liStlGJ tUSTAIiHfta» 15) 

CQMKjM/uTHck/ Mh, XVA*<»XUbWl 
GIMLfSiON MFU‘j,l‘^) t eFC:t»,l5 ) 

DO i: JJ=1.M 

Zu a . IIs:l#N 
I sIlH 

rf Ul, jJJsUd. j)-ALP 1* (lid ♦!, Jl * J»l - A».»’2»(U4it ♦UdtJ-lM 

♦►mlf.* alPx* lu ii*i, j-n*ur T^i * j* n>u(i»i, jd) ♦uM-itj-1) ) 

Ar (i^. t JJi - vd it 
23 rODllNUE 
1C rtNTMJE 

CALL XFUKiKNu fAFttiF) 
cu z: J=l ,N 

GO I=1,N 

UtMinroQ^T d./Ai.A>’DA(I, Jii 
VH=aF d.jj 
XHd, J)=DtM««* VH 
CoNfiMUE 
X.iEA'J- ij. 
nw j-i»M 

DO 5 ' 1 = 1 .N 
Xi-ltH i=XHLflN *XH d t Jl 
i>0 COMTIMUF 

X.-ic-A )?Xi:tAN/(N*NI 
TtrtP.sHAlC^FAlj 
T- Mf*i = PM:i*PA ‘i 

ai/A-«= ( l.-TrHPl)* d .-Tc.MP 2) ♦vAk/ U , *T £MPl » / d . ♦TEMP?) 


xC£\/i=SCkl IXVAki 
n £ T U 

£MH 


original i 




V-4 


SUftRQJll.^; QUANT 

r 0 n K U U / G 0 U ^ / U 1 1 7 ♦- 1 7 1 I P C 0 N t K U , , P A i U ♦ P A 0 1 1 A t PI f A LP2 , 
V ( 1 5 1 151 ,A*.A*’UA (lb, 15) , tfkP » tfil5 , 151 ,XH (15, 15), 
•XhSlly ,15) ,UuTAP( 15,15) ,uSTAP*hil5, 15) 

COHMtJN/CTHF P/ X.icAM , Xw A» , X 3E V 1 
tNAX = X ilt.ANt2.5*XGf VI 
Fi1lN = XMEAN-2.5*Xi)PVX 
OJ 3: J-ltN 
ao 3: i«t,N , 

NiJ=NLtV(I,J) 

ir(.«IJ ,un, C) GO TO G1 

M = ...... 

0= iFT.lAX-t.MlN) Z*^ 

. Vi=c..’l lN+0/2.,. . 

LX=»XHfi, JI -tltlNI/ O 
«HClO-V1UX*G 

iPlVtiJI^Ji .Ci,. FMAX) VHOLOsEMAX-0/2, 

IMXiiI,J) ,l£, LMINI V)iOLO=Vl 

GO n 50 

G1 VhuLJ -XMEAN 
•■O AnSii., JlstfHGLD 
53 CONI iNUt 
PLtU^N 
tND - .. 




iUtJiinUTiNc xFGRMlUfftOUM.BDUM J 

ulrtfc ’lSiuN tUUMliSriSI tUDUM ll5,li>l ( DJiUM 115 1 1& ) » 
•ZLIW (64) ,.'t (3l ,S(ai 

CG = -T^T 
IFSET = 1 

f,. ( 21 ;: r 

M 131-1 

cn 5 lST=ii2 
uO 5 1=1, N 

cn 2C . 1 = 1 , M 

JP = IM 1-1 1*^2 

' , , ■ ■ . . 

/ TjiJM ( 

ZUU'1 < jIls-0* 

IF (J ,Gl. 2 .AMu. J ,CF* N*il GC TO !«? 
un T*' 2C 

19 IFUiT .EO. ^1 A&u=AnuH(I, J-11 
lF(i5T .to. £) Ar,u = CDU«tI,j-l) 

2DUM(jP)=AbG 
22 rJMflNUE 

CALL HAFN(2DuM,M,INV,5,IFS£T,lF£KP 1 
err ?l J=1,MT 
jr = lMJ-l.»*2 

IFIJ .UC-* 2 .AUu. J .Lfc, N+l) GO TO 18 
bJ TO 21 
IF 

1FH-.T ,lP. ’» CPuHi J-i ,I>=MUO 
iFlI'=;r fetO. 21 oDU.1(J-l,i)-AvL) 

21 CONflri'J*' 

5 rjNTlMUt 
hetupm 

FNO 












SjWrOUTINi. fccon 

Cl> NHOm/GR nuh/u tl 7 , 1 7 > , H C Ct< , U 1 1 f N , P Al t t P A u 1 « aL P 1 »A LP2 


1&» »h,*.AKQAJ itfAR,Vti5,l5l ,XHU&*lS), 

♦XH4llt.,l'5> ,Ui,ThP. ilS.lS) ♦U5TARMI15, 15> 

COMMON/uT HcK/ KlitANtXVAKfXQF VI 

COHMv>H/OLiT/AMflN 

OIMEililON Af{15,16»,3Fti5,15) 

00. I"', Jsl»N,, ,, . . , . . 

CtMN-^QRTU,/A..M»*rjAII,jn 
UST4% tl, Ji-XH$ (I*J)*0£MiV 
AFUi JisUSTAnt tl.JJ 
10 CONllMUt. 

K.3-N 

call XFOfthiNC tAF »0F» 

DO J= 1 , h 
OG 3 : Irl,N 

USTA J)=BF a*4»>AK£AN 
30 CONTIVUF 
RTru^^'r : 








hisfo ;: / y \ y ; ; 

n 1 HP i S ION It! t btf ) , TdNP ( 6b » r S f MROL i 1 5C » » IG KAF M 1 1 

COMilU.t/ANAL YS/lJ(25bl T 

OaTh 1CT/6‘»‘* - / 1 TbliP/bS^a . / »ST-R/ tp* ni.«N</l H / 

•/ ^ ■■ V ■ 

iM'?GSR*2 lU . V "■ 

ktAL *G X y ' ■ V ^ ^ 

FOK'n T C////Z) I 

FU‘?«AT tiHi V ■ ■ I 

F jSMhT IIX, »rtAXlMaM = t,i.l2,5, lXrtrtlN|trtUM“t»£12,5, 
♦X, *TT?.a-»«hi2»5,i.X,»rtrAi'l = *,Fl2.5*iA, tSTA^4u».RD 
Fui^il&l (IXttPGuUUlit ,5X»>L0WLk»»5X,»QCUUitR^ 

♦c. StlNTii tt 16 f ♦ aCtUkFE.I4CF5*/iX,>LE(/iut t6X* tLtVfeU 
FokMAT t JtX,l£) i4X,l6l ) 

furmhI ax,ij?Ai) I 

FuRMfiT (uX, I2,tX,Ei3.b,iXt Ib,lX, ♦vt,iCCAl) 

FuRMxT (1 X t tt-ic-AN roUAkc tRFOR-trhlZ .|t23X, fROOl 

FOkMATClX, rSiGNML TO NOiSc RxTiO 1=|, ?1? i 5,17X# 
-no 2 = i-,El2.‘^.» : . - yyl ^ 

FOkSAmx.tSlGNAL TO MUITE PATIO I jlN OG SCAtE= 
*TU f<Oi^E I' AMO ? IN Lii ^ii,kLi.=>»ri2 
Format (IX, AO EOUaRE cRkLR XI Ll'3 SLALc = 

♦UA^t tRROR ll‘i. ut j CALt = ♦ » c 12 • a » , . [ 

x»iAx=:. . , f 

< .‘I T^'~ ^‘5b-» I 

iu‘n = j . 

' , ■"■ ■ . . ■ ‘ \y: 

r.RRSJ 1 = C. 

DO 1 1=1 » 255 
F5AD (A ) X 
RtAD( :<) lU 
DO 17 K=l ,25b 
;.;X(K) - lUlK) 
rONTiNUc 
DO 2 J=i, ?55 
aUMl= SUMl f X|J 0 
SJ42=0UK2*XU) *X 101 

tkRSOOi LftRSlJKt (A tjJ>-XXiO» 1 MX ( J) -X XUn 
IF t<i J) .uT.XfAXi XMAa=X(OI 
IF (X( J) .Lt.XOIOJ y«IN = X I J» 

LGUTIMUE ' ■ 
rONTIMUt 
K.AiW£= 127. 

:',T.nTftl-SUWl , . ■ ■ 

XXi1dAO=iUMl/l255.-255. > 


C» *0ARM32) *XXIII56) 
, DASH/ IH -/ , P„U^/ IH* 


IX, tRANG £“ » , Eli • 5 1 1 
D£VlAf 10 N=>, £12.51 
21X,>bHCH STAR REPR 


.It AN SQUARE tRRORs* 
♦SIGN At TO NOISE RA 
♦,€12.5, 5X, ♦SIGNAL 
5,UX,^ROOT MEAN SC 






£;=;>«l=“RRSUK/{i55,-?»r,, ) 

~ K ? * «/ 0 (t T ( E KR 1 ) 

SNl--?ftNGE/E;AR2 
5.Nl=$a/E,.f':i?>:; :: 
t tl :rli: .^Al OGa 
t t«'io , ♦ALCOi: : £,xH2) 
bNaHt-2C.^fiLOGlf: (SMi ) 
bNDii’’ = >G .•fi^noK; (SN2) 

101 ’ f^ANCic .TOTAii X Xf<Ffrf!,i,| 

KkI^^r icp, rkf?i,:FK«? 

iGi 

FkIkT ;illV EPktDci , EivRijB2 : ■ 

J- P INT 101 

•^klNJ r iC<5t b^ll»S^3 . :.^ 
pRt'jr t(:i 

1:10., :SNU:bi'^ 

REkPiD u 

XTW^ = lXKAX-XivlN)/Lt«, 

TtMP ( i ) =XHIN 
00 i ,i -=1,255 
kcAU {!♦ J X 

OO J=l,255: , 

00 3 <=1,04 

r£.M(- ( < f i) i TthF ( K'l ♦ X IMC 
IE (XCj|-l = N!P(Kfi)» 0,0,5 
.rtlNTTsjlUfc. 
icri'<<=icr(K)»i 
CCNTI^iUt 
rGmri;iut 
Ii.'rk,X'CO' ■■■■■' 
uo 7- 1=1,04,.:.;. 

r uilif 

NQsIJTMflX/io'] .♦!, 

PRI 4 T 101 

1C4, No' ^ 

.on ^ , 1=1 ,:iri ■:::;.. ''- y ' , ; : 

TGRA*';! { I ) =I'*’Mu*l '■ 

CO.MTINUE ; ^ - -. 0 

F=>iMr :.iQi ^ ^ 

P^^INT 1C5, iGFAPh 
no i'l 1=1,1 jr 





14 CCMTIMLir 

DO lb Xssrf isCiie 
Mktli -i'uUS 
to UJNTIMUE 

H‘?iNl tCo, OAF 

FNOsMJ 

Dft b I=l»f>4 

i Flirt (U .bE.NO) GO in g 
■Do: i: -Ksl vltC- 
!5YMjJL IK) sRLAr!K 
to rONIlMUE 
GO TO tt 

9 JrICTUI/PNU 

DO It K-t.J 
SYM3ou (K) sblAH 
12 COMi ’■'lUE 


on 1 ". K=j4.ft«t: 

SYMjj^iM sauAriK 
t& CuNTUUL 

tt >rinr iv7. I.itMF ti» tlDWI) , iSVf^OQLCKI 1,10 0 

0 CONTINUE 

CRlNT iSb. QAft 
F^1^4r ICF. IbKAPH 
!?cTU>N 


